



service_effect_sensitive = function(.df, .subset = NULL, ...){
  
  if (.subset== "all"){estimation_sample = .df}
  if (.subset== "java"){estimation_sample = .df %>% filter(java_indicator == 1)}
  if (.subset== "nonjava"){estimation_sample = .df %>% filter(java_indicator == 0)}
  
  bws <- seq(0.25, 10, by = .25)
  for(i in 1:length(bws)){
    
    temp_sample =
      estimation_sample %>%
      filter(abs(forcing_final) < bws[i])
    
    temp_sample$zscore = 
      kling_index_t1(.data = temp_sample,
                     .impute = F,
                     ... = ...)
    tmp = 
      temp_sample %>%
      lm_robust(zscore ~ pass_final, data = .) %>%
      tidy() %>%
      mutate(bw = bws[i])
    
    if(i == 1){fin = tmp} else{fin = bind_rows(fin,tmp)}
  }
  
  fin = fin %>% filter(term == "pass_final")
  return(fin)
}

java_index <- service_effect_sensitive(.df = data_4, .subset = "java", ... = q17_survey_java1, q17_survey_java2)
nonjava_index = service_effect_sensitive(.df = data_4, .subset = "nonjava", ... = q17_survey_java1, q17_survey_java2)
regional_index = service_effect_sensitive(.df = data_4, .subset = "all", ... = q18_survey_daerah1, q18_survey_daerah2, q18_survey_daerah3)
religious_index = service_effect_sensitive(.df = data_4, .subset = "all", ... = q19_survey_relg1, q19_survey_relg2, q19_survey_relg3)
national_index = service_effect_sensitive(.df = data_4, .subset = "all", ... = q20_survey_pancasila_h3, q21_survey_nation_ethnic_h3)
corruption_index = service_effect_sensitive(.df = data_4, .subset = "all", ... = q16_survey_transparent_h2, q13_survey_success_reason_1_merit_h2, 
                                            q13_survey_success_reason_3_sara, q13_survey_success_reason_2_connection_h2,
                                            q14_survey_test_vs_connection)




make_service_sens_plot <- function(df){
  
  df %>%
    ggplot(aes(x=bw, y = estimate)) +
    geom_point() +
    geom_line(color = "red") +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, color = "black") +
    ylim(-0.6, 0.6) +
    xlab("Bandwidth") +
    theme_bw() +
    geom_hline(yintercept = 0, linetype = "dotted") +
    scale_x_continuous(breaks=c(2, 4, 6, 8, 10),
                       labels=c("2","4","6", "8", "10")) +
    ylab("Estimate (SDs)") +
    theme(panel.grid.minor = element_blank(), 
          legend.position = "none",
          #axis.line.x = element_blank(),
          axis.ticks.x = element_blank(),
          panel.grid.major.x = element_blank(),
          #axis.line.y.left = element_blank(),
          axis.line = element_line(colour = "black"),
          panel.border = element_blank()) %>%
    return()
  
}

java_index <- make_service_sens_plot(java_index)
nonjava_index <- make_service_sens_plot(nonjava_index)
regional_index <- make_service_sens_plot(regional_index)
religious_index <- make_service_sens_plot(religious_index)
national_index <- make_service_sens_plot(national_index)
corruption_index <- make_service_sens_plot(corruption_index)


ggsave(filename = "./_4_outputs/figures/figure_a9a.png", plot = java_index, width = 6, height = 4)
ggsave(filename = "./_4_outputs/figures/figure_a9b.png", plot = nonjava_index, width = 6, height = 4)
ggsave(filename = "./_4_outputs/figures/figure_a9c.png", plot = regional_index, width = 6, height = 4)
ggsave(filename = "./_4_outputs/figures/figure_a9d.png", plot = religious_index, width = 6, height = 4)
ggsave(filename = "./_4_outputs/figures/figure_a9e.png", plot = national_index, width = 6, height = 4)
ggsave(filename = "./_4_outputs/figures/figure_a9f.png", plot = corruption_index, width = 6, height = 4)


